In the mean field description of superconductors in a magnetic field due to Abrikosov 
[1] , the zeros of the superconducting order parameter (vortices) form a triangular lattice in 
the low temperature ordered phase. This phase has a broken translational symmetry which 
reflects long range positional order, and phase coherence or off-diagonal long range order 
(ODLRO). The conventional theoretical picture [2] supposes that the two-dimensional 
Abrikosov vortex lattice will undergo a continuous phase transition at a melting temper- 
ature Tm, due to Kosterlitz-Thouless (KT) unbinding of paired dislocations in the 2D 
crystal [3]. However, it has been argued [4] that ODLRO in the Abrikosov state is de- 
stroyed by phase fluctuations due to thermal excitations of the lattice. Assuming that the 
loss of ODLRO is accompanied by the loss of positional order over the same correlation 
length scale £, the implication then is that only the disordered, vortex liquid state, will 
exist above zero temperature [5]. 

There have been various experimental papers which claim to support the KT melting 
scenario for thin films [6-8] . These all used different techniques, with changes in a measured 
quantity, such as resistivity or ac penetration depth, near the expected phase boundary 
given as evidence for the melting transition. The resistivity which falls rapidly towards 
zero as the temperature is lowered below a certain point is usually interpreted as being 
due to the formation of a flux lattice. However, this rapid change can also be explained 
as a consequence of increasingly long relaxation times which follow an Arrhenius law as 
the temperature is lowered [7] . Furthermore, in certain regions of the H-T phase diagram 
for BSCCO crystals, the high anisotropy make the separate layers behave as decoupled 
2D systems. Here the resistivity drop can also be explained by an activated process [9] 
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without any need to invoke a phase transition. 

Previous Monte Carlo (MC) studies of the 2D vortex system [10-13] have produced 
conflicting results. Tesanovic and Xing [10] find evidence for a transition, Kato and Na- 
gaosa [11] and Hu and MacDonald [12] conclude that the vortex lattice melts through a 
first order phase transition (for which there is no experimental support whatsoever) , while 
O'Neill and Moore [13] find no evidence of a vortex lattice at nonzero temperatures . 

This study extends the results of [13] by encompassing much larger systems and longer 
MC sampling times (made possible by a much improved algorithm), but nevertheless comes 
to the same conclusions, i.e. the absence of any phase transition at finite temperatures. In 
addition we can relate the relaxation times which we find in the system at low temperatures 
to the behavior of the resistivity as seen in Ref. [9] . 

Our starting point is the Landau-Ginzburg free energy functional of the complex order 
parameter ip(r), which can be written as [13] 

FW(r)]/k B T e = J rf 3 r{«(T)|^| 2 + /3(T)|^| 4 /2+|D^| 2 /2^}, (1) 

where D = —ihV — 2eA, the mean field zero-field transition temperature is T c , and ot{T), 
(3, fi, are the usual phenomenological parameters. 

The functional integral for the partition function can be approximated by expand- 
ing the order parameter in a basis set consisting of only the lowest Landau level (LLL) 
eigenstates <p p , of the gauge invariant operator D 2 /2//, i.e. 

N 

ip = Q^2v p (p p , (2) 

where Q = (3>o/ f3Bd) 1 ^, (with d the film thickness, B the external magnetic field), and 
N is the number of vortices [13]. The vortices are constrained to move on the surface of 
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a sphere of radius R which is imagined to contain a monopole of fixed strength at the 
center, such that 4nR 2 B = iV$ . The flux quantum $ is related to the magnetic length 
l m via $o = 27r/^S. Henceforth, Z m , which is the characteristic distance between vortices, 
is used as the measure of length; for example in these units, R = (N/2) 1 / 2 . An advantage 
of using this spherical geometry is that the vortex system has full rotational invariance - 
a feature absent in the geometry used in Refs. [11,12]. 

We neglect fluctuations in the A field and choose the gauge and eigenfunctions <p p to 
be the same as those in Ref. [13]. Within this LLL approximation, the free energy in Eq. 
(1), becomes 

N 2N 

F[{v p }]/k B T c = a T M 2 + iN- 1 \U P \ 2 , (3) 

where 

N 

U p = Y, h p . q h q B^\N - q + 1, q + l)0(p - q) v p _ q v q (4) 

q=0 

with the effective temperature variable «t = dQ 2 (a + eBK/n), 0(p — q) is the Heaviside 
step function, h p = [(N + 2)lpl/2nN(N -p)l} 1/2 , and B(x,y) = (x-l)!(y-l)!/(x + y-l)! 
is the Beta function. Large negative values of «r correspond to low temperatures , and 
mean field behavior is recovered (at 'zero temperature') in the limit «t — >• — oo. The LLL 
approximation is expected to be valid at large fields approaching the upper critical field 
H C 2, which correspond to small values of «t- 

The partition function requires integration of exp{— F[{f p }]//c B T c } over the set of 
complex coefficients {v p }. The form of Eqs.(3) and (4) permits an efficient MC sampling 
procedure to be devised. 

To examine vortex correlations, we evaluate a 'structure factor' defined on the sphere 
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as 

sr = ^J dndn' (mn)\ 2 \m')\ 2 )Yr(wrm, (5) 

where O is the solid angle, Y" z m are spherical harmonics, and < .. > denotes thermal 
averaging. A simpler expression can be calculated by assuming that in thermal equilibrium, 
the structure factor is azimuthally invariant, so that we may choose S(l) = Sf : and this is 
given by 

V 

where I p \ = Nh 2 j dQ \(p p (Q) | 2 P^(cos 9) exp (/ 2 /4), and Pi is the Legendre polynomial. 

The standard Metropolis algorithm has been used to sample the distribution of the 
dynamical variables {i> p }, with their initial configuration generated randomly. Typically, 
the first ten percent of the total number of Monte Carlo steps (MCS) used in each run 
are discarded to allow the system to equilibrate, and similarly, measurements from which 
averages are produced, are taken every ten MCS. System sizes of N = 120, 160, 200, and 
400 were studied. 

Some physical quantities such as the entropy (proportional to iV -1 J2 p {\v p \ 2 )) reach 
apparent equilibrium very rapidly (on the order of 5 x 10 3 MCS), and differs only by about 
one percent from the largest system (N = 400) to the smallest (N = 120). In comparison, 
much longer sampling times are needed to accurately determine the structure factor S(l) 
of Eq. (6), ranging from 8 x 10 5 to 2 x 10 6 MCS. In particular, for low temperatures 
(«t < —6), at least 10 6 MCS are required to obtain reasonably stable data. We will 
return to a more detailed discussion of the relaxation time scales later. 

First, we discuss the possibility of a first order transition. Numerical simulations in 
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the case of 2D systems have always produced conflicting results when searching for KT 
crystal melting [3]. Conclusions are often heavily dependent on factors such as (dislocation) 
core energies and are highly sensitive to boundary conditions. It is always difficult to do 
simulations for sufficiently large systems to minimize boundary and finite size effects. 

Following Refs. [11-12], we have searched for signs of hysteresis behavior in the energy 
and entropy, but without success. It is useful to investigate the probability distribution 
of the energy P(E), when attempting to interpret the nature of a phase transition [15]. 
Fig. la is a plot of this distribution, which shows no sign of the 'double peak' expected for a 
first order transition at the transition temperature reported in [11] and [12], corresponding 
to ax — —10, for up to N = 400. Only a single peak is observed (at all temperatures 
surrounding olt = —10) which narrows with system size at a fixed temperature . In the 
asymptotic limit, it is expected that the width at half height should tend to zero like 
iV -1 / 2 as N — > oo, and this is seen in Fig. lb. Both of the studies [11] and [12] employ 
quasi-periodic boundary conditions on the plane, such that the ground state is constructed 
to be a lattice. A possible reason then, for the appearance of a double peak, is that as the 
correlation length £ over which crystal order exists becomes of the order of the system size 
L, as it will for sufficiently low temperature, their system may easily fall into an apparent 
low energy state. However, at fixed temperature , this situation might change as L tends to 
infinity, although the crossover might only happen at very large system sizes. Furthermore, 
there is no sign of a first order transition in any experimental data that we are aware of. 

The structure factor of Eq. (6) as the temperature is lowered acquires peaks at values 
of I corresponding to the reciprocal lattice vectors G = |G| of the triangular lattice. For a 
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triangular lattice, the first peak is at G ~ 2.694, in terms of the normalised 'momentum' 
l/R. The scaling argument stemming from a zero temperature transition [5,13] predicts 
that the peak heights in the structure factor of the disordered or liquid phase, should 
scale as S(G) ~ £ 4 ~ |o;t| 4 - We have determined S(G) by extrapolating between two 
measurements of S(l) taken at l/R values on either side of I = G. That £ varies as \cct\ as 
a T — > — oo, was established in Ref. [13], and this has been reinforced by our own results 
[5]. In particular, the peak widths, which vary as l/£, decrease as |q;t| _1 as ar — > — oo. 

Turning now to the dynamics, the running value of the structure factor S(l) at the 
reciprocal lattice, / = G, was investigated systematically as a function of time, and we 
shall call this function S(G, t). This particular value of / was investigated in most detail as 
for this wavevector, the 'critical slowing down' associated with the zero temperature phase 
transition, is most pronounced and hence its study determines the longest relaxation time 
for the system. We find that sometimes, there appears to be more than one characteristic 
time scale. This is manifest by a local (in time) minimum in the value of S(G, t), which is 
different to the eventual true equilibrium value. However, the initially relaxed local value 
of S(G,t) only changes by a few percent as it 'creeps' to the next locally stable value, as 
seen in Fig. 2a. Any local minima clearly correspond to metastable states of the vortex 
system, the first of which might be reached after a relatively short time (< 5 x 10 5 MCS). 
In cases where even the initial relaxation time is long (> 5 x 10 5 MCS), or when S(G,t) 
slowly tends to another metastable state, we believe that the underlying mechanism for the 
longer timescales involves creating and moving a dislocation through correlated crystalline 
regions of the vortex system. The energy of such a dislocation [3] within a 'lattice' of size 
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£ is given by E dis ~ c 66 rf/^ log(£// m ), where c 66 is the shear modulus. Hence the activated 
relaxation time r ~ exp(£ , t ji s //csT) can be written as 

t = T exp{A|a T | 2 log(C|« T |)}, (7) 

where To, A, and C are constants, and we have used the relations \ar\ ~ £, ot\ ~ 
c 66 dl 2 m /k B T c [13]. 

The value of r has been extracted from the MC data by fitting to plots of S(G, t) 
at times after any initial relaxation, either growing ~ (1 — e _t / r ) or decaying ~ e _t / r 
exponentials tending to an estimated (possibly local), equilibrium value for S(G,t). This 
procedure is intended to measure the long timescale so that in cases where S(G, t) drifts 
away from the initially relaxed state at say time Tj to the next local minimum over a 
fitted time of r', we take the measured relaxation time to be r = Tj + r'. A more formal 
determination of r through a calculation of the autocorrelation function, was not feasible 
given the large relaxation times involved. 

The plot in Fig. 2b shows that Eq. (7) fits the data for a T < —4. The points for 
iV = 120 begin to bend below the straight line at very low temperatures, or large \ar\, 
and this is probably due to finite size effects (as can be seen from the trend for the larger 
systems), when the correlation length £ becomes comparable with the linear system size 
R. 

Yazdani et al [8], argued that by examining ac resistivity at high frequencies it is 
possible to probe superconducting films on length scales which are not affected by disorder. 
The samples used were MoGe films, and the large decades of response frequency covered 
(over a small temperature range) imply a spread of relaxation times much wider than 
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that obtained from our MC simulations. This is because disorder does dominate in the 
low frequency regime (corresponding to length scales larger than the Larkin-Ovchinikov 
length), at least in these conventional superconducting films. In the case of BSCCO 
crystals, where disorder is less important, high fields applied perpendicular to the Cu- 
O layers make the layers decouple and hence act as 2D systems over a portion of the H-T 
phase diagram [16], providing a more relevant comparison with our results. Resistivity 
measurements at various fields as a function of temperature below T c , show a gradual 
drop, typically exhibiting a 'knee', as the temperature is lowered. Below this 'knee', the 
resistivity can be explained as a field dependent activation process [9]. The relaxation 
process of Eq. (7) which is the so called 'plastic' timescale [16], is of this same Arrhenius 
form. Taking the resistivity to be inversely proportional to the relaxation time of Eq. (7), 
i.e. p ~ 1/t, it can be shown that our results match qualitatively the data of Briceno et 
al [9] over a temperature range from olt — —5 to —8, corresponding to the data points 
just below the knee. At lower temperatures, disorder is expected to be important again, 
leading to processes like TAFF which have timescales that dominate over the 'plastic' ones 
[16]. However, we have not been able to obtain a quantitatively accurate fit as it appears 
that the resistivity data in Briceno et al does not belong to the 2D scaling regime; in 
other words the measurements at different fields do not collapse onto a single curve as 
a function of ax, implying that the data is outside the LLL description. (This is not 
surprising, given for instance, the small region in the H-T diagram which fit the 2D scaling 
regime in magnetisation data for BSCCO [17]). It would be useful therefore, if a systematic 
investigation of the resistivity was made in the regime where the layers are decoupled and 
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where the LLL approximation is valid. This would provide a direct comparison with the 
simulation results reported here. 

We would like to thank Y. Kato, N. Nagaosa, Jun Hu, A. Yazdani, A.D. Rutenberg 
and M.A. Pettifer for helpful discussions. The numerical simulations were performed with 
DEC Alpha 3000/400 workstations and on a Cray EL98. This work was supported through 
a SERC grant. 



Figure captions 

FIG. 1. (a) The energy distribution curve, P(E) plotted against E/Emf where 
Emf = — N\ar\ 2 /2/3a is the mean field energy (with /3a — 1-159 the Abrikosov fac- 
tor), at «t = —10 for N = 400, showing only a single peak. The lack of a double peak 
suggests the absence of a first order phase transition, (b) The peak width at half height W, 
as a function of A -1 / 2 . The straight line extrapolates to zero, indicating the asymptotic 
regime has been reached. 

FIG. 2. (a) Running value of S(G, t) as a function of MC run time, for N = 200 
at «t = —8 showing an initial relaxation to a local metastable state after 8 x 10 5 MCS, 
and then the subsequent 'creep' to another metastable state after 1.3 x 10 6 MCS in to- 
tal, (b) The logarithm of the relaxation time for various system sizes plotted against 
|o;r| 2 log(|aiT|). A linear fit is obtained for temperatures lower than olt — —5, verifying 
Eq. (7). The falloff at larger \ar \ are due to finite size effects. 
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